Entangling many-body bound states with propagative modes in Bose-Hubbard systems 
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The quantum evolution of a cloud of bosons initially localized on part of a one dimensional optical lattice and 
suddenly subjected to a linear ramp is studied, realizing a quantum analog of the "Galileo ramp" experiment. 
The main remarkable effects of this realistic setup are revealed using analytical and numerical methods. Only 
part of the particles are ejected for a high enough ramp, while the others remain self-trapped. Then, the trapped 
density profile displays rich dynamics with Josephson-like oscillations around a plateau. This setup, by coupling 
bound states to propagative modes, creates two diverging condensates for which the entanglement is computed 
and related to the equilibrium one. Further, we address the role of integrability on the entanglement and on the 
damping and thermalization of simple observables. 
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The last decade breakthrough experiments on ultra-cold 
atomic gases have revived the field of strongly correlated 
many -body quantum systems UJ, especially within non- 
equilibrium aspects Q. The very low dissipation rate and 
long time phase coherence of these systems allow the inves- 
tigation of genuine quantum effects, the role of integrability 
on thermalization properties, and the possibility to engineer 
desired states. Two examples are particularly relevant for this 
study. The first is Bloch oscillations (BO) |3| which occur 
when a particle travels on a lattice and experiences a con- 
stant external force F (potential ramp): Its momentum drifts 
with time according io q{t) = q{0) -\- Ft (we set h and the 
lattice spacing equal to one), but modulo the Brillouin zone, 
setting the BO period tb = 27t/\F\. They have been ob- 
served in many physical domains : semiconductors (H, ther- 
mal gases f5l, photonics |6|, cold atoms |7 1, and phonons f^. 
BO can survive to the many -body regime, with a damping 
possibly related to the integrability of the model |[9l[T0 |. The 
second example is the release of atoms from a trap, a stan- 
dard protocol with cold atoms. Yet, keeping the optical lat- 
tice on during the expansion allows one to handle metastable 
states 1 1 1 1, or study transport phenomena in disordered poten- 
tials 1 12 |. We lastly stress that ultra-cold atoms, beyond their 
fundamental interest, could be used as devices to mimic optics 
experiments |13 | or semi-conductor physics |14|. 

In this paper, we put forward an experimentally realistic 
setup which ejects interacting bosons living on an optical 
lattice using a linear potential (see Fig. [T]). This "Galileo 
ramp" experiment displays remarkable features that are un- 
derstood with analytical and numerical calculations. Only 
part of the particles is ejected and forms a wave-packet which 
shape and number of particles are determined. Thanks to the 
initial correlations, these traveling particles remain strongly 
entangled with the ones remaining self-trapped in the initial 
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region, hence creating two diverging and entangled many- 
body condensates. In the self-trapping region, particles ex- 
hibit Josephson-like oscillations reminiscent of BO, with a 
density plateau due to many -body interferences, and damp- 
ing controlled by the integrability of the model. The setup 
is particularly versatile, in comparison to analog proposals to 
create wave-packets |T51, and relies on the general idea of 
coupling bound states to the propagative modes of a lattice. 
Results in the hard-core bosons limit provides a quantitative 
and intuitive understanding which allows for straightforward 
generalizations (mirrors, beam-splitters,. . . ). 



I. SETUP DESCRIPTION 

A. Hamiltonian 

We describe the setup by means of the one-dimensional 
(ID) Bose-Hubbard model (BHM) 

3 3 3 

where h^- is the bosonic creation operator at site j, rij = h^-hj 
is the density operator. J and U are respectively the hopping 
and interaction magnitudes while Vj (t) is the time-dependent 
external potential. At times t < 0, the potential is a deep box 
of width A which confines all N particles in region A, with 




FIG. 1. (color online) Sketch and notations of the considered setup. 
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the average density p = N/A. At times t > (quench), the 
potential is suddenly changed to a linear ramp Vj =Va — Fj 
in region A and zero elsewhere (see Fig. [T]). We provide in 
Appendix |A| the criteria to neglect tunneling toward the upper 
Bloch band. 

For simulations, the setup is embedded in a larger box of 
size L, and the gas is assumed to be isolated. Its dynamics 
is studied via analytical methods in the non-interacting and 
Hard-Core Boson (HCB) U = oo limits. HCB dynamics 
can be solved exactly using the mapping to a fermion Hamil- 
tonian and diagonalizing the single-particle physics (see be- 
low). Observables are then expressed as determinants that 
are computed numerically. The out-of-equilibrium calcula- 
tions are based on Refs.[T6land[T7l Otherwise, two ab-initio 
techniques are used to compute observables: Lanczos diag- 
onalization and time-dependent density-matrix renormaliza- 
tion group (tDMRG) ifTsHIH (see also Refs. 22 for a good 
introduction), which is particularly useful to compute the en- 
tanglement entropy of Fig. |4] Lanczos diagonalization are 
used with a cutoff M in the onsite boson number. M = N 
is taken on the box configuration (except for L = 18 for 
which M = 7). For L = 32, the 64bits limitation im- 
poses M = 3 < N. This can participate in the little dif- 
ference between the free bosons analytics and the numerics 
in Fig. igd). tDMRG is performed with M = 2. After 
sweeping through the chain to apply the evolution operator, 
we keep the maximal value of the block entanglement entropy 
Smax {t) and we update the number of kept states m at time 
t -\- dt to m(t + dt) = mt[m{t) x ex.p ASmax{t)], in order 
to account for the eventual growth of entanglement. Typi- 
cally, starting the evolution with m(0) G [16,24], we finish 
with m(T) e [24, 110] depending on the value of U and F. 
Comparisons with Lanczos calculations ensured that this is 
sufficient to monitor errors for the considered parameters. For 
these ab-initio methods, the Hilbert sizes at play are not chal- 
lenging and the numerical errors are under control. 



movers, to which are attributed the local velocities v{±q{x)), 
where q{x) = arccos[(y(x) — e)/A]. This approach allows 
us to reconstruct the evolving density profile p{x; t). 



II. WAVE-PACKET EMISSION OF THE GALILEO RAMP 

A. Wave-packet density profile 

The typical evolution of the density of HCB after the 
quench for increasing forces is given on Fig.|2ja). Only part 
of the bosons are ejected, with maximum velocity A, and the 
wave-packet spreads with time. In the HCB hydrodynamic 
description, the population with energy £ G [—A, A] is con- 
nected to the propagating states outside A (see Fig. [2jb)). 
These particles escape from region A and propagate toward 
oo. As for BO, the momentum of a right(left) mover of 
energy e obeys q^{t) = ±q -\- Ft. From energy conser- 
vation, it converts all its potential energy into kinetic en- 
ergy when it reaches the position A at a time given 
by Ft^ = zfq -\- arccos(— e/A). After escaping from A, 
both travels at constant velocity, giving the ballistic trajectory 
x^{t]q,e) = VA^ - s'^{t-t^)-^A. At times larger than r^, 
these particles have all left region A and their density profile 
is the sum of the right and left movers contributions 

P^c(^;^)= / ^ / ^^{x-x^{t',q,€)) , (1) 
J -A F Jq^s) 27r 

in which Q{€) refers to the yellow domain of Fig. [2jb). 
Fig. [2jb) shows the very good agreement, up to small quan- 
tum interference effects, between exact diagonalization data 
and the hydrodynamic prediction ([T]). 

B. Number of ejected particles 



B. Hydrodynamic description of HCB 

In the limit N^A ^ 1, after a Jordan- Wigner transforma- 
tion, a chain of HCB maps to the free fermion Hamiltonian 
V, = ^q^qVlVq with the fermionic operators r]q of momen- 
tum q G [— TT, tt] . The kinetic part of the energy, — A cos q with 
A = 2 J, builds up a ID band of width 2 A and gives a velocity 
v{q) = A sinq. The ground-state is the corresponding Fermi 
sea of Fermi momentum qr = 7rp. This is a good description 
of the initial state provided the well is large and deep enough. 
The features of the density dynamics are well captured by the 
hydrodynamic limit using the continuous position variable x. 
The initial state has approximately a coarse grained phase- 
space density ^(x, ^) = ^n[o,A](^)n[_g^,g^]((7), where 
Ii[a,b] is the characteristic function over the interval [a, 6]. Af- 
ter the sudden quench, the energy of a particle at position x is 
shifted by the additional potential value V{x). The dynamics 
being unitary, quasi-particles are emitted to the right and to 
the left on a trajectory of constant energy e. The density of 
particles is the sum of the densities of these right and left 



The total number N^^^ of escaping HCB is readily ob- 
tained by integrating the initial density on the same domain 

p{x < A,e;t = 0) = [ — 5{€ - V{x) ^ Acosq) . 

One obtains that, for a small positive force F < A(l + 
cos qF)/A,8ill particles are connected with propagative states, 
i.e. N^^^ = N. When increasing F beyond, some particles 
remain trapped in region A and the explicit calculation yields: 

^^^"^ = ^ [-/F + Sin 9f) - 9{FA/A)] (2) 

with g{x) — x(2 — x) — {x — 1) arccos(x — 1) when 
A(l + cos(7f)M < F < 2A/A, and g{x) = when 
F > 2 A/A. A comparison with simulations is displayed on 
Fig.[2jc), demonstrating the accuracy of the approach. 

Lowering interactions naturally modifies A^esc- Taking the 
opposite case of free bosons, this single-particle physics is 
solved along the same lines : free bosons are initially gathered 
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FIG. 2. (color online) (a) time-evolution of the density profile at half-filling for increasing ramp height Va = 1, 3, 5, 8 [J], (b) sketch 
explaining the hydrodynamic calculations. Inset: numerical escaping profiles (full lines) compared to the hydrodynamic approach (dashed 
lines), (c) number of escaping HCB particles vs. Va for various initial densities p (hydrodynamic predictions in full lines), (d) same number, 
varying the interaction strengths [// J at half-filling, and compared to analytical limiting cases. 



in the ground-state wave-function of the deep well so that the 
initial density profile in region < x < Ais 
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In this regime, all bosons share the same dynamics and the 
density profiles of the trapped and outgoing wave-packets de- 
pend only trivially on N which enters as a prefactor. The 
initial state thus consists of N bosons of initial momentum 
qo tt/A. The particles leaving the confinement region are 
again those with energy e G [—A, A], corresponding to the 
space interval i3 = [max(0, A — 2A/F), A]. Then, one ob- 
tains N^lf = NforO<F< 2A/A and 
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for F > 2 A/A. This time, the strong force scaling is F~^. 
This last result is also valid for A/" = 1, in which case the 
single-particle wave-function is in a superposition of bound 
and diffusive states. In Fig. [2jd), numerical simulations for 
various U shows how results interpolate between these two 
limiting cases with an interaction-dependent behavior at large 
F. Note that, surprinsingly, interactions do not favor escaping 
at small forces. 



III. DYNAMICS OF TRAPPED PARTICLES 

A. Density profile 

As seen from Fig. [2ja), the trapped HCB exhibit a rather 
rich physics with oscillations reminiscent of BO. This can be 
understood in the hydrodynamic approach: We consider the 
main contribution coming from trapped particles with ener- 
gies e G [A, Va — A] (grey area in Fig.[2];b)). Neglecting 
tunneling escapes at the band edges, the corresponding densi- 
ties of right/left movers are 
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V{x) ^ Acos{q±Ft)) 



These are clearly periodic functions of time, oscillating with 
period tb- By introducing q = q ± Ft in each term and inte- 
grating over energies, the coherent part of the density reads 

Pcoh{x;t) = / ^[AyA-A]{V{x) - Acosq) . (5) 

JFt-qp 

From ([5]), we see that oscillations are present provided p <1 
(superfluid regime of HCB). Indeed, at unit filling p = 1 (Mott 
state) qr = and the integration range then covers a full pe- 
riod, leading to a static trapped density profile. The density 
profile can actually be computed at any time t from ^ (see 
below), and a comparison with numerics is given in Fig.jSja). 
We first observe that incoherent contributions due to high en- 
ergy modes of energy e G [Va — A, Va + A] (green area in 
Fig.|2jb)) not taken into account in ^ increase both the aver- 
age and the fluctuations of the particle density on the left side. 
Indeed, due to the free boundary at x = 0, high energy parti- 
cles sharing initially the same momentum q are not reflected 
at X = at the same time. Consequently, a dephasing appears 
between them and their contribution to the total spatial den- 
sity is somehow incoherent. A simple way to suppress this 
incoherent effect is to eject these particles into a propagative 
band on the left (x < 0), as shown in Fig.jSjb) where we see 
a very good agreement, up to small interference effects, be- 
tween the exact numerical results and the hydrodynamic pre- 
dictions when the incoherent particles are removed. 

Another striking feature of Fig. |3ja)-(b) is the existence of 
a stationary density plateau: The argument of the IT func- 
tion in ([5]) reaches its edges for positions XM/sup given by 
V{xm) — Acosq = Va — a and V(xsup) — Acqs ^ = A. 
For maxxinf < x < minxsup, the integral is independent of 
time and the initial density plateau survives. Therefore, the 
condition to have this plateau is F > 4 A /A and its width 
is A — 4A/F. The boundaries of the plateau are given by 



Xp = max Xinf and x^ 



min Xc„n, which can be written 
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At each sides of the plateau, an excess density oscillates, being 
at the left(right) for (half-)integer multiples of tb • 
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FIG. 3. (color online) Density profile without (a) and with (b) a propagating band opened at the left of x = at different times for Vj^j J — 10, 
A — 200 at half-filling p — 1/2 (hydrodynamic predictions are given by the dashed red lines), (c) Corresponding current density at x = A/ 2. 
Inset: Current profiles at different times (colors refer to the circles). Full lines are the exact results while dashed red lines are the hydrodynamic 
predictions, (d) same but for a smaller ramp Va/ J — The straight red line is the expected truncated sine function while the dashed red one 
is the pure sine oscillation (see the main text). 



The explicit density profile is easily computed at any time 
t from Eq. ([5]). For example, at integer multiples of the BO 



period r^, the distribution of the trapped bosons is shifted 
maximally to the left (down to Xinf(g'F) < x~) and is given 
by 



Pcoh{x]nTB) 



\ arccos[(V(x) - Va)/A 4- 
P 

p - ^ arccos[V(x)/A - 1] 




1] < X < Xir^Mr) 

X+ < X < Xsup(gF) 
XsupiQF) <X<A 



(6) 



The distribution then propagates from this macroscopic left 
state to the macroscopic right one (reached at half periods 
times t = {n -\- 1/2)tb with integer n). The spatial profile 
of the right state is simply deduced from the left one by the 
transformation x A — x, using the mirror symmetry at 
A/ 2. These are the results plotted in Fig.jsja-b). 



B. Collective Josephson-like oscillations 

This remarkable pendulum motion of the density is natu- 
rally associated with a flow of particles through the plateau 
region. The flow of particles from the left to the right of the 
trapping zone and vice versa gives rise to a periodic current 
density j{x;t) which, in the hydrodynamic limit, is simply 
the sum over all momenta of p{q)v{q), i.e. the quasi-particle 
currents [p'^{q) — p~ {q)]Asmq. In the trapped region the 
current density is given by the integral expression 

Ft+qp 

j{x;t) = A J ^ smqUiA,VA-A]{V{x)-Acosq). (7) 

Ft-QF 

For F > 4 A /A, in the plateau region Q one has simply from 
^ a spatially constant current (which was expected from the 
continuity equation since the local density is constant in time) 

j{x; t) = — sm{qF) sin(Ft) , (8) 



which oscillates harmonically with period tb • The maximum 
current amplitude is obtained at half filling when qp = 7r/2, 
while the current is exactly zero in the Mott phase (unitary 
filling giving qp = tt) reflecting the stationarity in time of the 
trapped density in that case. These Josephson-like oscillations 
are reported in Fig. [3jc) and we see that the hydrodynamic 
description matches perfectly the numerical data. Outside the 
plateau, the current density has a spatial variation which is 
simply deduced from its integral representation ([7]). Again, as 
seen in the inset of Fig.jSjc), the hydrodynamic description is 
very good. 

At lower forces, for 2 A/ A < F < 4:A/A, when there is no 
more a plateau region but still self-trapped particles, the tem- 
poral evolution of the current is no more given by a pure sine 
function due to the interplay between the integration interval 
and the support of the door function entering into For ex- 
ample, at half filling qp = 7r/2 and taking Va = 3 A, one has 
from ^ in the middle of the ramp 



f sin(Ft) for 
A/27r for 



sm{Ft)\ < (Va/2A- 1) 
sm{Ft)\ > (Va/2A- 1) 



This truncated sine essentially reflects the escape of particles 
which were initially located in the middle of the condensate. 
Indeed, for Va < 4A the escape locus x+ extends over the 
middle of the ramp, since x+ < A/ 2, leading to a lowering 
of the current intensity after the corresponding particles have 
been ejected. As a support of this, notice on Fig.jSjd) how dur- 
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ing the first half period, as the particles are moving from the 
left to the right and have not yet abandoned the condensate, the 
current shows a perfect sinusoidal signal which gets truncated 
as time goes on. Notice that, despite a behavior similar to the 
Josephson effect, the oscillations stem from a fundamentally 
different effect - many -body interferences with strong interac- 
tions - and without a tunneling barrier - the plateau develops 
within the system. 



IV. ENTANGLEMENT ENTROPY BETWEEN THE 
WAVE-PACKETS 

After the quench, the condensate is split into two entan- 
gled pieces moving apart: the escaping particles and the 
self-trapped ones. This entanglement between bound states 
and propagative ones can be quantified through the bipartite 
von-Neuman entropy St{A) = —TY{gt{A)\ngt{A)}, where 
gt{A) = Tr>A{|^(t))(^(t)|}. As seen on Fig.|4];a), this en- 
tropy essentially evolves up to an asymptotic value 5*00 {A), as 
expected, which depends on N, U and F. The effect of the 
force is qualitatively inferred from Fig.[4ja). At small F when 
all particles are ejected (see F = J/ 4:), after an initial increase 
of the entanglement due to the crossing of position x = A 
by the many -body wave-packet, the entanglement finally van- 
ishes when most particles have left region A. This small-F 
regime is almost independent on the interaction. Increasing 
the force, the asymptotic entanglement first increases, passing 
through a maximum when approximately half of the initial 
density is ejected, before decreasing at larger forces F, simply 
because A^esc falls down. The effect of interactions on Soo {A) 
is two-fold : one is to modify A^esc as seen on Fig.[2jd), and the 
other one is related to the fact that non-integrability (signifi- 
cant when U J) usually increases the chaoticity of excited 
states, leading to higher entanglement entropies. A strong en- 
hancement of the entanglement is thus observed for U = 2J, 
by comparison to the behavior close to the integrable HCB 
limit (Fig. |4ja)). Lastly, in the inset of Fig. |4ja), we show 
the early evolution of St{A) which is controlled by the rate 
at which particles are emitted. This initial emitting rate, as 
expected, is enhanced by the repulsion strength U and by the 
force F. 

The strong entanglement observed between the two pieces 
of the condensate must stem from the initial state in which 
all particles, embedded in the same condensate, are naturally 
entangled. Therefore, the asymptotic value (A) of the bi- 
partite entropy should be related to the entanglement present 
initially between the particles that will leave region A and 
those that will stay. In the HCB limit, this idea actually gives 
a quantitative prediction for Soo{A). Indeed, in Fig.Qb) we 
show 5*^=800 (^) — ^oo(^), as a function of Va, compared to 
the initial equilibrium ground- state entropy So{x^ (Va)) eval- 
uated at the point x+ (Va) (see Fig.|2];b)) which marks the left- 
most initial position of the ejected particles. The entropy 
is computed from exact diagonalization and also plotted from 
the (continuous limit) conformal field prediction 1231 
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FIG. 4. (color online) (a) evolution of the entanglement entropy 
St{A) for an increasing force F {p — 0.5, A = 8, L = 32; full lines 
U / J — 2, dashed lines U / J — 20; Inset: short time evolution.) (b) 
Asymptotic value of St {A) (taken at t = 800) for HCB as a function 
of Va and compared with the initial bipartite entropy Sq{xp {Va)) 
from both numerics and CFT predictions (p = 0.5, A = 200). 



with d ^ 0.25. The good agreement found between these 
three curves supports the above picture. The remarkable con- 
nection between the equilibrium and non-equilibrium entan- 
glements found in this setup could be helpful in the search 
for measuring many -body entanglement This connection 
could be possibly used to infer equilibrium (local) correla- 
tions from measurements on far apart particle packets. We 
also stress that the setup works for few-particles physics for 
which the measurement of entanglement is easier. 



V. PROPOSALS ON ACHIEVING WAVE-PACKET 
MANIPULATIONS 

In Fig. [5] we give ideas on how to manipulate travelling 
wave-packets generated by the previously discussed "Galileo 
ramp" setup. The arguments are intuitively based on the HCB 
picture from energy conservation and demonstrated by nu- 
merical calculations. First, a travelling wave-packet can be 
catched on a lattice by suddenly increasing/decreasing the 
chemical potential (in a ramp shape on the figure) so that the 
particles cannot escape the trapping zone. It requires that the 
condensate is on the region before changing the potential (first 
sketch). If the potential is raised before the condensate arrives, 
then the region acts as a mirror (second sketch). If the mirror 
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height is finite and a barrier builds up, we obtain a situation 
similar to single-particle tunneling process. Then, the barrier 
acts as a beam-splitter (third sketch). Last, applying a non- 
uniform potential on a wave-packet can reshape it by slow- 
ering the fastest particles and accelerating the slowest ones 
(fourth sketch). This would help fight the natural broadening 
of the wave-packet. In order to show that these intuitive be- 
havior does take place in the real time evolution of HCB, we 
give in Fig.[5]two examples of numerical simulations showing 
the realization of a mirror and of a beam-splitter. These de- 
vices, together with the proposed source of entangled wave- 
packets, could help realize interferometry measurements and 
their applications - with the advantage that the two wave- 
packets can be spatially well separated, thus being able to ex- 
perience different real- space paths. 



VI. DAMPING AND THERMALIZATION IN THE CLOSED 
BOX GEOMETRY 

We see that the main physics is well captured by the in- 
tegrable HCB and non-interacting limits. Still, the model is 
essentially non-integrable at finite U/ J, which has some con- 
sequences on the dynamics. In order to focus on this aspect, 
we remove the propagative band and confine the particles in a 
box, thus allowing larger sizes for simulations. In this config- 
uration, BO are damped in the chaotic regime |9 1 and experi- 
ments reported very low damping nearby the non-interacting 
point ifTOl . Here, we go further by computing the damping 
time and the distance from thermalization for all U. In particu- 
lar, we stress that the integrable nature of the HCB limit shows 
up in both quantities. Therefore, we relate in this section the 
spectral features of the final Hamiltonian and of the quench 
distribution to experimentally accessible observables. In par- 
ticular, we would like to see whether the (non) -integrable na- 
ture of Hamiltonian has noticeable consequences on simple 
observables. In addition to the damping of BO which was pre- 
viously discussed in Refs.[10|and|24| we address the question 
of the thermalization |i2l- 



A. Time evolution after a sudden quench 

The system is prepared in the ground- state |?/^o) of a box 
without the linear potential which is suddenly turned on at 
t = 0. We write the energies and \n) the eigenvectors of 
the final Hamiltonian governing the dynamics. In this basis, 
the matrix elements of an observable O are Onm = {m\0\n) 
and the associated excitation frequencies between n and m 
are written ujnm = En — Em- Let = (^|V^o) be the co- 
efficients of the initial state in this basis, and Pn = |cnP the 
corresponding weights. These weights are the quench distri- 
bution, or diagonal ensemble, and simply provide the aver- 
aged contribution of excited states to the time-evolution. The 
fidelity J^, defined by 

Ht? = \A{t)\^ with A{t) = {mm = ^Pne-*""«* , 

n 

plays an important role in the dynamics and contains the infor- 
mation on the distribution through the Lehmann representa- 
tion of A{uj). Given an Hermitian observable O, its real-time 
behavior and corresponding derivative read 

0{t) = O + ^ 2\CnC''^Onm\cOs{uJnmt + 0nm) (10) 
n<m 

0{t) = - ^ 2\CnC^0nm\^nm Sm{uJnmt + (/>nm) (H) 
n<m 

with 

= Y,PnOnn (12) 

n 

the time-averaged expectation value (assuming a non- 
degenerate spectrum) and (/)nm some phases determined by 
the initial conditions and the observable. There are two in- 
teresting behaviors in the time-dependence : the short-time 
behavior, related to the damping of the observable towards O, 
and the value O itself which can be compared to statistical 
ensemble predictions to check thermalization. While only the 
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FIG. 6. (color online) (a) Evolution of the momentum distribution nq(t) as a function of time (vertical shift) for different U/J and a fixed 
force F = 2 J. The damping seems to be clearly related to the non-integrability of the model, (b) normalized damping of the momentum 
distribution ng(t) for different sizes and interactions (the signal local maximum is taken over a short time window to remove finite size 
oscillations associated with the ^-space discretization). 



diagonal elements enter the expression of O, the off-diagonal 
ones Onm will play a role in the damping, together with the 
associated weights Ic^c^l and frequencies cj^m- These two 
features should depend on the integrability of the Hamiltonian 
through selection rules, for instance. 



B. Observables 

In this section, we define the observables O that are used. 
The simplest observable is the local density nj{t) = {nj){t). 
Computing the one-body density-matrix gjk{t) = {h^jhk){t) 
allows one to access two important observables. The first one 
is the momentum distribution nq{t) = {nq){t), with 

«.W = ^Ee^'^'"'WW (13) 

jk 

which is measured using time of flight techniques. A peak in 
Uq usually measures the coherence of the condensate at the 
corresponding wave- vector. The second one is the condensate 
fraction fo{t), defined as the largest eigenvalue of the one- 
body density-matrix. It matches N when all bosons share the 
same state, and is of order ^/N in the HCB regime. Even so 
it is not experimentally observable, the fidelity T{t), defined 
above, quantifies the distance from the initial state and is as 
well sensitive to the chaotic features of the Hamiltonian. 



C. Damping of observables 

The spectral features and chaoticity of the Hamiltonian are 
provided in Appendix [B] When one goes to the large- /7 or 
large-F limits, the frequencies spectrum Unm is such that 
they are all very close to the same level spacing (U or F), 
i.e. a nearly equally spaced spectrum. Of course, as judged 



by the poissonian nature of the spectrum at finite J, a dense 
level- spacing distribution is still present, generically leading 
to damping. In this regime, we expect some observables to 
show strong oscillatory behavior at the main frequency and 
its harmonics. This is a situation qualitatively related to the 
classical motion in closed orbits. In the quantum version, 
this would show up in the pseudo phase-space (O, O) (from 
Eqs. ([TO|-([TT])) as more or less complex orbits depending 
on the number and weights of the harmonics involved. In 
the nearly equal level spacing spectrum situation, the fidelity 
would obviously display similar resonances (or revivals). In 
the case of non-integrable Hamiltonian, although level repul- 
sion tends to rigidify the spectrum, we expect that the weights 
and the frequencies are spread and of the same magnitude 
so that the damping of observables usually occurs on shorter 
times. In the pseudo phase-space (0,0), damping translates 
into a spiral structure of the orbit. We do observe this qualita- 
tive effect, particularly on the local density and its associated 
particle current (data not shown). 

The first interesting quantity to look at is the momentum 
distribution v.s. time because it captures well the Bloch os- 
cillations. We recall that, at the single-particle level, the mo- 
mentum is shifted with time according to q{t) = q{0) -\- Ft. 
In Fig.[6j we recover that nq{t) is globally shifted with time 
for all U/J, with a period close to the Bloch period tb- In- 
creasing U further does not change much this period (hardly 
seeable in the numerics) and the main effect of interactions 
is to destroy the main peak signaling the bosonic coherence. 
Indeed, for U = 2 J, we observe a strong damping of the co- 
herence after short times. Interestingly, increasing U brings us 
close to the HCB limit and the damping time increases again. 
Einally, for very large U, the central peak shape is preserved 
for very long times. 

Another interesting point with the momentum distribution 
within this setup is that the zero-momentum evolution is qual- 
itatively related to the fidelity of the system in the near- 
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FIG. 7. (color online). Comparison between riq^o (t) and the fidelity 
in the integrable regimes, showing the revivals. 




10 20 

t [n/j] 

FIG. 8. (color online). Typical evolution of the condensate fraction 
fo{t) with time for various interactions and fixed ramp F = 2J. 
Inset: inverse damping time r v.s. U/J. 



integrable regimes, both signaling strong revivals of the initial 
wave-function (see Fig. [7]). In fact, we know that the mo- 
mentum distribution gets back to its initial position in g-space 
after a period tb . More than the momentum distribution, the 
whole many-body wave-function actually comes back close 
to the initial state, yielding strong revivals in the fidelity. We 
observe that the peaks in the fidelity and in n^^o are nicely 
correlated signals in these two limits. 

Fig. [8] shows the damping of the condensate fraction /o for 
increasing U/J and fixed force. The behavior is very similar 
to the momentum distribution, but easier to fit to extract a typi- 
cal damping time r. We extract the typical damping time from 
the normalized condensate fraction time-evolution by fitting 
the curves using the function /(t) = a + (l — a)/(l + (t/r)^). 
a is the asymptotic value of the fonction, and r is such that the 
curve has decayed of half its distance from 1 to a. The expo- 
nent b is found to vary between 1 and 5 in the fits but its value 
does not significantly change r which is the quantity of inter- 



est. Good fits are obtained at small U while oscillations spoil 
the fit at large U, which nevertheless provides a reasonable 
estimate of the damping time. The main plot and the inset 
supports that the rs become much longer near the two inte- 
grable points. 



D. Thermalization 

In addition to the short-time behavior, the difference be- 
tween integrable and non-integrable regimes can show up in 
the long-time average value O. We see that both the fea- 
tures of the diagonal distribution pn and the behavior of the 
observables with energy intervene in Eq. ([12]). If one wants 
to study thermodynamics only, the energy distribution pn 
is the only useful quantity. In the case of a generic dis- 
tribution, the equivalence of ensemble should be sufficient 
to provide the same thermodynamics as usual thermal en- 
sembles 1251 . If one wants to compare time-averaged ob- 
servables O to those obtained from thermal ensembles, the 
hypothesis that O nn ^ ^(-^) ha.rdly varies within the en- 
ergy shell given by the typical energy fluctuations naturally 
leads to the identification with O and, therefore, to ther- 
malized observables in this sense. Justification of this idea 
(sometimes dubbed as the Eigenstate thermalization hypoth- 
esis, or ETH) has been proposed from quantum chaos argu- 
ments i26l|23 and studied numerically |28, 29 1. Some ex- 
tensions and more details prescription of this statement have 
been discussed more recently 130 IK stressing the importance 
of the correlations between pn and Onn- In this context, 
integrability or the proximity to integrable points on finite 
systems has been shown to lead, most of the time, to non- 
thermalized observables QT] |28l l30ti35l . The particular sit- 
uation of an interaction quench in the Bose-Hubbard model 
which possesses integrable limits, has been discussed in depth 
in Refs. l3Ql[32ir36 , and J7, Of course, the above approaches 
suffer from limitations in their applicability. For instance 
it is expected to work best at high energies where quantum 
chaos describes well the states. Yet, this corresponds to high- 
temperature where the physics is in general not so interesting. 
Another issue comes when taking the thermodynamic limit 
N ^ oo after or before time-averaging. 

We now turn to the calculation of the time-averaged expec- 
tation O of an observable compared to its thermal prediction 
O*. O is directly obtained from ([12]) from full diagonaliza- 
tion of the Hamiltonian on small systems (we use L = 12 
and = 6) and from computing the quench distribution p^ 
exactly. is obtained using calculation in the canonical en- 
semble. A microcanonical description is also possible but it 
is in practice less transparent as one has to tune to energy 
window. Usually, the outcome does not depend much on the 
choice of the distribution at sufficiently high energies (tem- 
peratures) and looking at simple observables; Besides, the 
equivalence of ensemble suggests that both ensemble should 
in principle yield the same results in the thermodynamical 
limit, i.e. for large enough systems. The temperature Tb of 
the Boltzmann distribution is directly obtained by demanding 
that the mean-energy should be the same as the mean-energy 
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FIG. 9. (color online). Comparison between time- averaged profiles obtained after the quench and thermalized profiles with the same mean- 
energy (in a closed box). Left: for a force F large enough to put enough energy in the system, but not too large to destroy non-integrability, the 
gas looks thermalized in the non-integrable regime U / J — 2.0, while clear deviations from thermalization are observed close to the integrable 
limits U — 0, oo. Right: When the force is larger, this finite-size system {L — 12) looks too integrable to observe approximately thermalized 
behaviors. In both cases, Tb indicates the temperature obtained from the mean-energy using the canonical ensemble. 



of the system after the quench {E) 
Boltzmann weights read Pn 
tion function and we plug them in ([12]) to get O*. 



Y^nPnEn- Then, the 
/Z with Z the parti- 



The goal is now to see whether deviations from thermaliza- 
tion are related to the integrability of the model studied pre- 
viously. In Fig. [9j we show the profiles for the three simplest 
observables for quenches with moderate forces F and three 
typical values ofU/J respectively close to the non-interacting 
limit, deep in the non-integrable region and close to the HCB 
limit. One must not choose a force which is too small, other- 
wise little energy is put in the system, the temperature would 
be too small and the discussion possibly spoiled by finite-size 
effects |32|. If the force is too large, we know from Fig. 12 
that the system looks integrable for all U/ J. We find that 
F/Jc^0.5,lisa reasonable compromise to use the quench 
protocol to probe the relation between thermalization and in- 
tegrability in this setup. Indeed, in Fig.[9j we observe that for 
F/J = 0.67, all profiles agree well with the thermal predic- 
tions for [//J = 2, deep in the non-integrable regime. Close 
to the integrable points, deviations are found but the effect de- 
pends on the chosen observable. We find that the local kinetic 
energy gjj-\-i is the most suitable observable to probe the ef- 
fect of the integrability (something also discussed in Ref.l33]). 
For F/J = 1.67, it is much harder to reach thermalization, 
probably due to the near integrable features of the Hamilto- 
nian for this small system. 




100 



FIG. 10. (color online). Relative distance (in %) between the time- 
average and the thermal expectations of an observable v.s. U/J (den- 
sity rij, kinetic energy Qjj+i, momentum distribution Uq). f mea- 
sures the non-integrability of the model. 



relative distance from thermal prediction and dub it: 
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(14) 



In order to quantify the deviation from the thermal predic- 
tion, we average over the sites j (or momentum q for n,) the 



This is the quantity plotted in percentage in Fig. 10 as a func- 
tion of U/J. We qualitatively expect that non-thermalized 
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regimes occur close to integrable points fSS", 31 -'331. This is 



well observed in Fig. 10 for which the integrability measure, 
denoted by f (see Appendix |B] for definitions), is manifestly 
correlated to the non-thermalization measure. The advantage 
of this setup with respect to interaction quenches 1 3^ ^32, .36| 
[37l is that the ratio U/J is kept fixed. In other words, it de- 
couples the quenching parameter from the parameter which 
mainly controls integrability. 



VII. CONCLUSION 

We introduced and characterized a remarkable setup, the 
"quantum Galileo ramp", which can be used to create en- 
tangled many-body wave-packets. We further provide sim- 
ple views on how to manipulate the condensate afterwards. 
An important result is that the asymptotic non-equilibrium 
entanglement entropy is quantitatively related to the equilib- 
rium initial one in the HCB limit. In addition, non-trivial 
Josephson-like oscillations are found and we show the setup 
is well suited to study the role of integrability on damping and 
thermalization. We underline that results on HCB concerning 
the local density are also valid for free fermions. Beyond the 
field of cold atoms, the setup could be developed for polariton 
condensates 1,38 J . 



Acknowledgements 

This work is supported by ANR-09-BLAN-0098-01. M. C. 
benefited from the International Graduate College on Statisti- 
cal Physics and Complex Systems between the universities of 
Nancy and Leipzig. 



Appendix B: Spectral features and chaos in the box geometry 
(without a propagative band) 

1. Hamiltonian 

We study using Lanczos and full diagonalization techniques 
the Bose-Hubbard Hamiltonian in box of size L and subjected 
to a linear potential Vj = Va — Fj, where F is the slope (or 
force): 



n 



j 



-E 

j 



(Bl) 

with Vj the external potential. If the onsite boson cutoff M = 
1, we have the integrable XX model while if M = N we 
recover the Bose-Hubbard model. In the following, all data 
are for a density p = N/ L = 1/2 (superfluid regime) and full 
diagonalization are done with L = 12 (Hilbert space size of 
about 12000 states), and Lanczos calculations up to L = 18. 

This model has been investigated in a similar context, ana- 
lyzing the effect of interactions on BO and its regular/chaotic s 
regimes. The main known results are that a new period, in ad- 
dition to tb, emerges in the BO for finite U and in the strong 
force limit, which reads 2it/U |9|. Chaotic motion is found 
for small F and in the presence of interactions, leading to a 
damping of the BO 1241 . A strong reduction of damping has 
been observed experimentally fTOl when ^ 0. The Mott 
regime U ^ J of the BO has also been briefly investigated 
in Ref . 40 Level statistics and the chaotic nature of the spec- 
trum has been partially analyzed in Ref. 24 and 41 Approx- 
imate methods such as discrete non-linear Schrodinger equa- 
tion and mean-field theory have been also used to study the 
dynamics of this model |i40l |42l . We lastly mention that 
the hard-core boson dynamics in a tilted bichromatic lattice 
(which possesses a gap in the single-particle dispersion rela- 
tion) has been investigated in Ref. ,44| 



2. Density of states 



Appendix A: Neglecting the tunneling toward the upper Bloch 

band 



Experimentally, tunneling to upper Bloch band could occur 
if the band is to close in energy to the states of the trapped par- 
ticles. Indeed, under the constant force F, Bloch oscillations 
may be ruled out by Landau-Zener interband transitions (and 
particles may possibly escape into the continuum if the upper 
gaps are too small). The rate at which this escape toward the 
upper band is expected is roughly given by the Landau-Zener 
tunneling formula e~^^ 1^ |[5] [39l where ^ is the interband 
energy gap and c is a constant depending in particular on the 
recoil energy {F^i = A here). To avoid the escape of parti- 
cles toward the upper band one needs accordingly to apply a 
sufficiently low force such that F <C cb^ . On the other hand, 
one needs a sufficiently high force in order to scan the first 
Brillouin zone within a reasonable small time (from an exper- 
imental point of view) 



We first discuss the nature of the many-body density of 
states (DOS) of the Hamiltonian as a function of the two 
parameters \J j J and F j J . There are three energy scales in 
the problem : the hopping J which sets the bandwith of the 
single-particle energies, the local interaction IJ and the ex- 
ternal potential ramp F. These typical energies are visible 
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in the many-body DOS of the spectrum as seen in Fig 
When only J is present, the DOS has a gaussian shape cen 
tered around zero energy. In the large-F limit, we observe 
equally-spaced peaks separated by the energy F. This struc- 
ture is reminiscent of the single-particle Wannier-Stark ladder 
spectrum. It is as well trivially understood at the many-body 
level when F ^ J since F corresponds to the potential en- 
ergy cost when a particle jumps from one site to its neighbor. 
Similarly, when F = and U ^ J, elementary processes 
corresponding to changing the onsite number of particles yield 
an equally-spaced spectrum of energy U (Mott lobes). No- 



tice that F rapidly kills these lobes as we see on Fig. 11 that 
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FIG. 11. (color online). Evolution of the density of states of the 
many-body Hamiltonian with slope F for different U . We observe 
main bands separated by U at large [/ or by F at large F . 



F/ J > 1 is sufficient to destroy the lobes in the DOS. 



3. Signatures of integrability and chaoticity in the spectrum 

In order to study the non-integrable nature of the Hamil- 
tonian ( [Bl| ) which governs the time-evolution, we use level 
and wave-function statistics (not shown). However, we see 
that the usual unfolding procedure for the spectrum will be 
plagued in the large F or large U regimes, due to the peak 
structure of the DOS. Following Refs. 45 and 46, we use 
the statistics of the ratio of consecutive level spacings = 
min{5n,5n-i)/ max((5n, ^^n-i) where 5n ^n+i - is 
the level spacing and are the energies. Probability distri- 
butions P(r) are then compared to poissonian and gaussian 
orthogonal ensemble (GOE) statistics. 

An example for the model under study is given in Fig. [TT 
in which we observe both regular (Poisson) and chaotic 
(GOE) distributions, depending on the parameters. In or- 
der to extract a simple number to be plotted against the pa- 
rameters, we use the normalized mean ratio defined by f = 
((r) - (r)poisson)/((r)GOE - (r)poisson), which is 1 in the non- 
integrable regime and in the integrable regime. When F ^ 
0, the reflection symmetry of the box is lost and the Hamil- 
tonian ( |B1| ) has no spatial symmetries which could plague 
the level statistics if not taken into account. As discussed in 
Ref. 45 in the F = case, only the /7 = and J = lim- 
its are integrable, Bethe-ansatz predicting the Hamiltonian is 



non-integrable as soon as [/, J 7^ 0. Still, because of finite 
size effects, small and large U/ J regions look almost inte- 
grable on a finite cluster. The thermodynamic limit is par- 
ticularly hard to investigate numerically, although the trend 
is compatible with the Bethe-ansatz predicting the Hamilto- 
nian is non-integrable as soon as /7, J ^ 0. Still, because of 
(a) ^ 




(b) 1 



FIG. 12. (color online), (a) Probability distributions of the mean 
adjacent level spacing ratio P(r) as a function of U / J for two dif- 
ferent slopes F. (b) Evolution of the normalized mean adjacent level 
spacing ratio (r) as a function of U / J for different slope F. 



finite size effects, small and large U/ J regions look almost 
integrable on a finite cluster. The thermodynamic limit is par- 
ticularly hard to investigate numerically, although the trend is 
compatible with the Bethe-ansatz perspective (451. We then 
expect I24I im the large U and large F regime to be close to 
the integrable classical limit (J = 0), at least on a finite chain. 
The results for different F as a function ofU/J are displayed 

4 J, the distribution looks 



on Fig. 12 and we see that for F 
always poissonian. In order to observe non-integrable effects, 
we must not choose F too large (F < 1 for this p = 1/2 
and L = 12). We also notice from this figure that the optimal 
range of U to be in the non-integrable regime depends on F. 
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